network analysis raw
Introduction
Our philosophy in this iteration of the analysis is that we can simply treat the number of transects as a covariate to be controlled statistically. This avoids the need to bootstrap and throw away a lot of data, but we pay a price in terms of model identifiability due to the collinearity between elevation and transect number.
Prepare environment
Load data
Calculate network/group-level metrics
This is straightforward because none of these metrics involve independent abundance data (i.e. our survey data). We will focus on the network-level metrics of connectance and NODF (nestedness) and the group-level metrics of niche overlap and C score (both are measures of dietary overlap / partitioning)
Temperature ~ elevation
The question has been raised of whether we might learn more by using temperature instead of elevation as a predictor variable. We see here, though, that total degree-day accumulation and elevation are almost perfectly correlated, so I think there is little to be gained by using temperature as an explanatory variable. We do, however, see that 2011 and 2012 were substantially warmer than 2010. 
Extract species-level metrics
How well does mean proportional marginal total predict mean proportional flower cover?
It’s actually not so bad.
Since we have curvy relationships and nested variables, we will analyze our data using hierarchical generalized additive models (HGAMs)
I think there are two ways we could approach this.
If our goal is simply to explain network metrics, we can include all covariates (floral abundance, floral richness, BB abundance, BB richness, transect number) in the models, understanding that we have collinearity problems, and see which ones explain the most. Because of the (curvilinear) collinearity between elevation and the other covariates, though, this approach would make it harder to detect and interpret the infleunce of elevation, which is the main question we’re interested in.
The second approach is to say that while network metrics may, indeed, be influenced by abundance and richness covariates, these are just mediated effects of elevation. We may not care exactly how (i.e. via which mediating covariate) elevation influences network metrics; we might just want to know the shape of the response. But I think we must include transects as a covariate even in this approach; otherwise, we would need to go down the road of bootstrapping/rarefaction, and I think we all agree that we don’t want to do that.
My recommendation would be to go with option 2, and that is what will be presented below.
Network/group-level metrics
For these models, we will use the bs = “fs” factor smooth method. This preserves degrees of freedom by forcing all smooths to have the same wiggliness.
Results
- Niche overlap (15% explained)
- Linear negative effect of elevation
- No significant effect of transect number
- C score (30% explained)
- Positive quadratic/cubic effect of elevation
- No significant effect of transect number
Discussion These models have a lot of unexplained variation, but they suggest that the diets of co-existing bumble bee species become increasingly dissimilar (niche overlap) or disaggregated (C score) with increasing elevation. This pattern does not appear to be influenced significantly by transect count.
# Web asymmetry
gam_web.asym <- gam(web.asymmetry ~
s(elev.mean, year, bs = "fs", k = 5) +
s(bb.transects, year, bs = "fs", k = 4),
data = netmet,
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
check.gamViz(gam_web.asym) ##
## Method: REML Optimizer: outer newton
## full convergence after 14 iterations.
## Gradient range [-1.609747e-05,4.434315e-05]
## (score -42.69759 & scale 0.01290583).
## Hessian positive definite, eigenvalue range [1.090476e-06,36.6072].
## Model rank = 28 / 28
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(elev.mean,year) 1.50e+01 1.21e-04 0.74 <2e-16 ***
## s(bb.transects,year) 1.20e+01 5.52e+00 1.14 0.89
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(gam_web.asym)##
## Family: gaussian
## Link function: identity
##
## Formula:
## web.asymmetry ~ s(elev.mean, year, bs = "fs", k = 5) + s(bb.transects,
## year, bs = "fs", k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.50578 0.06483 -7.801 5.27e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(elev.mean,year) 0.0001214 14 0.00 0.544
## s(bb.transects,year) 5.5163473 10 10.58 7.83e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.592 Deviance explained = 62.3%
## -REML = -42.698 Scale est. = 0.012906 n = 74plot(sm(gam_web.asym, 1)) +
l_fitLine(alpha = 0.6)
plot(sm(gam_web.asym, 2)) +
l_fitLine(alpha = 0.6)
# Niche overlap
gam_niche.overlap <- gam(niche.overlap.HL ~
s(elev.mean, year, bs = "fs", k = 5) +
s(bb.transects, year, bs = "fs", k = 4),
data = netmet,
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
check.gamViz(gam_niche.overlap) ##
## Method: REML Optimizer: outer newton
## full convergence after 13 iterations.
## Gradient range [-4.387283e-06,8.27703e-06]
## (score -52.89884 & scale 0.01218919).
## Hessian positive definite, eigenvalue range [2.061461e-06,36.53777].
## Model rank = 28 / 28
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(elev.mean,year) 1.50e+01 2.32e+00 0.83 0.065 .
## s(bb.transects,year) 1.20e+01 2.56e-05 1.06 0.610
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(gam_niche.overlap)##
## Family: gaussian
## Link function: identity
##
## Formula:
## niche.overlap.HL ~ s(elev.mean, year, bs = "fs", k = 5) + s(bb.transects,
## year, bs = "fs", k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.28475 0.01284 22.18 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(elev.mean,year) 2.321e+00 14 0.73 0.00668 **
## s(bb.transects,year) 2.558e-05 11 0.00 0.86758
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.123 Deviance explained = 15.1%
## -REML = -52.899 Scale est. = 0.012189 n = 74plot(sm(gam_niche.overlap, 1)) +
l_fitLine(alpha = 0.6)
plot(sm(gam_niche.overlap, 2)) +
l_fitLine(alpha = 0.6)
# C score
gam_C.score <- gam(C.score.HL ~
s(elev.mean, year, bs = "fs", k = 5) +
s(bb.transects, year, bs = "fs", k = 4),
data = netmet,
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
check.gamViz(gam_C.score) ##
## Method: REML Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-5.627583e-06,3.739096e-06]
## (score -36.37735 & scale 0.01741642).
## Hessian positive definite, eigenvalue range [2.064089e-06,36.63284].
## Model rank = 28 / 28
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(elev.mean,year) 15.000 4.934 0.83 0.04 *
## s(bb.transects,year) 12.000 0.697 1.09 0.76
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(gam_C.score)##
## Family: gaussian
## Link function: identity
##
## Formula:
## C.score.HL ~ s(elev.mean, year, bs = "fs", k = 5) + s(bb.transects,
## year, bs = "fs", k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.43572 0.04316 10.1 4.09e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(elev.mean,year) 4.9343 14 1.379 0.000469 ***
## s(bb.transects,year) 0.6967 10 0.134 0.075160 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.238 Deviance explained = 29.7%
## -REML = -36.377 Scale est. = 0.017416 n = 74plot(sm(gam_C.score, 1)) +
l_fitLine(alpha = 0.6)
plot(sm(gam_C.score, 2)) +
l_fitLine(alpha = 0.6)
Species-level d’
First let’s look at the abundance and distribution of each species in each year. It’s probably not worth trying to fit the model to species that are extremely rare or limited in their distirbution. We will drop B. humilis and B. lapidarius from 2010, B. humilis, B. gerstaeckeri, and B. mucidus from 2011, and B. humilis, B. hypnorum, and B. mucidus from 2012.
bb_range.year <- net %>%
left_join(site_data) %>%
group_by(site, date, elev.mean, bb.sp, year) %>%
summarize(abund = n()) %>%
group_by(site, elev.mean, bb.sp, year) %>%
summarize(abund = mean(abund)) %>%
group_by(bb.sp) %>%
mutate(elev.floor = min(elev.mean),
elev.ceiling = max(elev.mean),
elev.range = elev.ceiling - elev.floor,
elev.med = elev.floor + ((elev.ceiling - elev.floor)/2))
ggplot(bb_range.year, aes(reorder(bb.sp, elev.med), elev.mean)) +
geom_line(size = 2, alpha = 0.25) +
geom_point(aes(size = abund), alpha = 0.5) +
facet_wrap(~year, ncol = 1) +
theme_light(12) +
scale_size_continuous(name = "Abundance") +
scale_color_discrete(name = "Abund. class") +
xlab("BB species") +
ylab("Elevation") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
For species-level d’, we will use the “by = factor” factor smooth method. This allows different wiggliness and independent significance tests for each species, which I think is the right way to handle species-level differences in the response of d’ to elevation. We will also create separate models for each year so that we can drop the rare species from each year and avoid a complicated interaction term (bb.sp * year) in the by=factor call.
Results (only significant terms shown)
- 2010 (33% explained)
- B. gerstaeckeri (positive, linear)
- B. pratorum (positive, quadratic)
- B. pyrenaeus (negative, quadratic)
- B. soroensis (positive, quadratic)
- 2011 (41% explained)
- B. hortorum (positive, linear)
- B. psithyrus (positive, linear)
- Transects (positive, quadratic)
- 2012 (37% explained)
- B. sensu-strictu (positive, cubic)
- B. hortorum (negative, linear)
- B. mendax (negative, leveling off above 1500)
- B. pascuorum (negative, linear)
- B. pyrenaeus (negative, linear)
Discussion The majority of species do not exhibit a significant relationship between elevation and discrimination, and which species do varies markedly by year. Only B. hortorum and B. pyrenaeus show significant relationships in more than a single year. Transect number was significant only in 2011. The most interesting finding is that the relationship between elevation and discrimination is generally positive in 2010 and 2011, but it changes to negative in 2012.
### Species-specific models
gam_d.2010 <- gam(d ~
bb.sp +
s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
s(bb.transects, k = 4),
data = filter(spmet, year == 2010 &
!bb.sp %in% c("humi", "lapi")),
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
check.gamViz(gam_d.2010) ##
## Method: REML Optimizer: outer newton
## full convergence after 16 iterations.
## Gradient range [-1.237157e-06,5.174296e-06]
## (score 31.42053 & scale 0.05353188).
## Hessian positive definite, eigenvalue range [4.464382e-07,83.51918].
## Model rank = 73 / 73
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(elev.mean):bb.spbss 4.00 1.00 0.92 0.140
## s(elev.mean):bb.spgers 4.00 1.00 0.92 0.145
## s(elev.mean):bb.sphort 4.00 1.00 0.92 0.120
## s(elev.mean):bb.sphypn 4.00 1.00 0.92 0.145
## s(elev.mean):bb.spjone 4.00 1.00 0.92 0.085 .
## s(elev.mean):bb.spmend 4.00 1.00 0.92 0.115
## s(elev.mean):bb.spmont 4.00 1.00 0.92 0.095 .
## s(elev.mean):bb.spmuci 4.00 1.36 0.92 0.150
## s(elev.mean):bb.sppasc 4.00 1.00 0.92 0.130
## s(elev.mean):bb.spprat 4.00 2.63 0.92 0.125
## s(elev.mean):bb.sppsit 4.00 1.00 0.92 0.135
## s(elev.mean):bb.sppyre 4.00 2.33 0.92 0.090 .
## s(elev.mean):bb.spsoro 4.00 2.35 0.92 0.115
## s(elev.mean):bb.spwurf 4.00 1.00 0.92 0.100 .
## s(bb.transects) 3.00 1.00 0.90 0.080 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
summary(gam_d.2010)##
## Family: gaussian
## Link function: identity
##
## Formula:
## d ~ bb.sp + s(elev.mean, by = bb.sp, bs = "tp", k = 5) + s(bb.transects,
## k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.35200 0.04632 7.599 2.25e-12 ***
## bb.spgers 0.05871 0.09037 0.650 0.5168
## bb.sphort -0.12442 0.07015 -1.774 0.0780 .
## bb.sphypn -0.01320 0.09957 -0.133 0.8947
## bb.spjone -0.01235 0.10967 -0.113 0.9105
## bb.spmend -0.24637 0.12576 -1.959 0.0518 .
## bb.spmont 0.08340 0.15245 0.547 0.5851
## bb.spmuci 0.08644 0.45429 0.190 0.8493
## bb.sppasc -0.10226 0.06966 -1.468 0.1440
## bb.spprat 0.05830 0.06713 0.868 0.3865
## bb.sppsit 0.15755 0.08034 1.961 0.0516 .
## bb.sppyre 0.01964 0.15497 0.127 0.8993
## bb.spsoro 0.09216 0.06646 1.387 0.1674
## bb.spwurf 0.08295 0.06551 1.266 0.2072
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(elev.mean):bb.spbss 1.000 1.000 1.297 0.25639
## s(elev.mean):bb.spgers 1.000 1.000 8.288 0.00452 **
## s(elev.mean):bb.sphort 1.000 1.000 1.784 0.18352
## s(elev.mean):bb.sphypn 1.000 1.000 0.697 0.40518
## s(elev.mean):bb.spjone 1.000 1.000 0.087 0.76805
## s(elev.mean):bb.spmend 1.000 1.000 3.378 0.06789 .
## s(elev.mean):bb.spmont 1.000 1.000 0.007 0.93355
## s(elev.mean):bb.spmuci 1.358 1.593 0.065 0.84897
## s(elev.mean):bb.sppasc 1.000 1.000 0.230 0.63212
## s(elev.mean):bb.spprat 2.632 3.157 3.208 0.02442 *
## s(elev.mean):bb.sppsit 1.000 1.000 2.882 0.09144 .
## s(elev.mean):bb.sppyre 2.326 2.693 3.655 0.02507 *
## s(elev.mean):bb.spsoro 2.345 2.849 2.965 0.04087 *
## s(elev.mean):bb.spwurf 1.000 1.000 0.808 0.36996
## s(bb.transects) 1.000 1.000 0.346 0.55741
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.196 Deviance explained = 33%
## -REML = 31.421 Scale est. = 0.053532 n = 196print(plot(gam_d.2010, select = c(2, 10, 12, 13)), pages = 1)
gam_d.2011 <- gam(d ~
bb.sp +
s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
s(bb.transects, k = 4),
data = filter(spmet, year == 2011 &
!bb.sp %in% c("humi", "gers", "muci")),
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
check.gamViz(gam_d.2011) ##
## Method: REML Optimizer: outer newton
## full convergence after 14 iterations.
## Gradient range [-3.815379e-06,1.975298e-05]
## (score -22.71315 & scale 0.03124357).
## Hessian positive definite, eigenvalue range [5.346636e-07,104.0106].
## Model rank = 68 / 68
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(elev.mean):bb.spbss 4.00 1.00 0.64 <2e-16 ***
## s(elev.mean):bb.sphort 4.00 1.24 0.64 <2e-16 ***
## s(elev.mean):bb.sphypn 4.00 1.34 0.64 <2e-16 ***
## s(elev.mean):bb.spjone 4.00 1.00 0.64 <2e-16 ***
## s(elev.mean):bb.splapi 4.00 1.47 0.64 <2e-16 ***
## s(elev.mean):bb.spmend 4.00 1.93 0.64 <2e-16 ***
## s(elev.mean):bb.spmont 4.00 1.00 0.64 <2e-16 ***
## s(elev.mean):bb.sppasc 4.00 1.24 0.64 <2e-16 ***
## s(elev.mean):bb.spprat 4.00 1.00 0.64 <2e-16 ***
## s(elev.mean):bb.sppsit 4.00 1.00 0.64 <2e-16 ***
## s(elev.mean):bb.sppyre 4.00 1.00 0.64 <2e-16 ***
## s(elev.mean):bb.spsoro 4.00 1.00 0.64 <2e-16 ***
## s(elev.mean):bb.spwurf 4.00 2.02 0.64 <2e-16 ***
## s(bb.transects) 3.00 2.44 0.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
summary(gam_d.2011)##
## Family: gaussian
## Link function: identity
##
## Formula:
## d ~ bb.sp + s(elev.mean, by = bb.sp, bs = "tp", k = 5) + s(bb.transects,
## k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.42663 0.03635 11.736 < 2e-16 ***
## bb.sphort -0.09794 0.05313 -1.844 0.06670 .
## bb.sphypn -0.25762 0.06042 -4.264 3.08e-05 ***
## bb.spjone -0.16543 0.05307 -3.117 0.00209 **
## bb.splapi -0.16576 0.06144 -2.698 0.00756 **
## bb.spmend -0.06793 0.09951 -0.683 0.49560
## bb.spmont -0.01154 0.05881 -0.196 0.84457
## bb.sppasc -0.02033 0.05378 -0.378 0.70575
## bb.spprat -0.08001 0.05182 -1.544 0.12417
## bb.sppsit -0.04304 0.05535 -0.777 0.43777
## bb.sppyre -0.01830 0.08372 -0.219 0.82719
## bb.spsoro 0.06489 0.05140 1.262 0.20829
## bb.spwurf 0.13937 0.05243 2.658 0.00847 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(elev.mean):bb.spbss 1.000 1.000 0.104 0.74756
## s(elev.mean):bb.sphort 1.241 1.441 7.304 0.00694 **
## s(elev.mean):bb.sphypn 1.339 1.599 1.047 0.43774
## s(elev.mean):bb.spjone 1.000 1.000 1.831 0.17750
## s(elev.mean):bb.splapi 1.472 1.789 1.856 0.24195
## s(elev.mean):bb.spmend 1.926 2.212 2.326 0.10547
## s(elev.mean):bb.spmont 1.000 1.000 1.828 0.17790
## s(elev.mean):bb.sppasc 1.236 1.429 3.808 0.05709 .
## s(elev.mean):bb.spprat 1.000 1.000 1.333 0.24958
## s(elev.mean):bb.sppsit 1.000 1.000 19.677 1.48e-05 ***
## s(elev.mean):bb.sppyre 1.000 1.000 2.695 0.10221
## s(elev.mean):bb.spsoro 1.000 1.000 2.279 0.13264
## s(elev.mean):bb.spwurf 2.021 2.490 1.427 0.21996
## s(bb.transects) 2.436 2.799 3.593 0.03613 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.319 Deviance explained = 40.8%
## -REML = -22.713 Scale est. = 0.031244 n = 235print(plot(gam_d.2011, select = c(2, 10, 14)), pages = 1)
gam_d.2012 <- gam(d ~
bb.sp +
s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
s(bb.transects, k = 4),
data = filter(spmet, year == 2012 &
!bb.sp %in% c("humi", "hypn", "muci")),
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
check.gamViz(gam_d.2012) ##
## Method: REML Optimizer: outer newton
## full convergence after 16 iterations.
## Gradient range [-1.188507e-06,8.313178e-06]
## (score -11.72191 & scale 0.03378739).
## Hessian positive definite, eigenvalue range [1.858343e-07,95.01981].
## Model rank = 68 / 68
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(elev.mean):bb.spbss 4.00 3.10 0.89 0.020 *
## s(elev.mean):bb.spgers 4.00 1.00 0.89 0.070 .
## s(elev.mean):bb.sphort 4.00 1.13 0.89 0.060 .
## s(elev.mean):bb.spjone 4.00 1.00 0.89 0.050 *
## s(elev.mean):bb.splapi 4.00 1.61 0.89 0.045 *
## s(elev.mean):bb.spmend 4.00 1.83 0.89 0.030 *
## s(elev.mean):bb.spmont 4.00 1.00 0.89 0.055 .
## s(elev.mean):bb.sppasc 4.00 1.00 0.89 0.035 *
## s(elev.mean):bb.spprat 4.00 1.54 0.89 0.035 *
## s(elev.mean):bb.sppsit 4.00 1.00 0.89 0.050 *
## s(elev.mean):bb.sppyre 4.00 1.00 0.89 0.020 *
## s(elev.mean):bb.spsoro 4.00 2.31 0.89 0.060 .
## s(elev.mean):bb.spwurf 4.00 1.00 0.89 0.055 .
## s(bb.transects) 3.00 1.00 0.86 0.020 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
summary(gam_d.2012)##
## Family: gaussian
## Link function: identity
##
## Formula:
## d ~ bb.sp + s(elev.mean, by = bb.sp, bs = "tp", k = 5) + s(bb.transects,
## k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.420864 0.036953 11.389 < 2e-16 ***
## bb.spgers 0.221879 0.096305 2.304 0.02234 *
## bb.sphort -0.085015 0.054129 -1.571 0.11799
## bb.spjone -0.211027 0.080103 -2.634 0.00914 **
## bb.splapi -0.159219 0.076893 -2.071 0.03978 *
## bb.spmend 0.210935 0.098327 2.145 0.03324 *
## bb.spmont -0.110492 0.070557 -1.566 0.11907
## bb.sppasc -0.024553 0.052687 -0.466 0.64175
## bb.spprat -0.077977 0.053391 -1.460 0.14585
## bb.sppsit 0.115933 0.056859 2.039 0.04288 *
## bb.sppyre 0.032636 0.088826 0.367 0.71373
## bb.spsoro 0.016465 0.052236 0.315 0.75296
## bb.spwurf -0.003091 0.052170 -0.059 0.95282
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(elev.mean):bb.spbss 3.098 3.591 4.855 0.00272 **
## s(elev.mean):bb.spgers 1.000 1.000 0.355 0.55181
## s(elev.mean):bb.sphort 1.134 1.254 5.310 0.01242 *
## s(elev.mean):bb.spjone 1.000 1.000 0.188 0.66530
## s(elev.mean):bb.splapi 1.614 1.961 0.344 0.68781
## s(elev.mean):bb.spmend 1.830 2.198 3.308 0.03324 *
## s(elev.mean):bb.spmont 1.000 1.000 0.491 0.48435
## s(elev.mean):bb.sppasc 1.000 1.000 4.411 0.03704 *
## s(elev.mean):bb.spprat 1.545 1.895 0.855 0.49405
## s(elev.mean):bb.sppsit 1.000 1.000 0.694 0.40587
## s(elev.mean):bb.sppyre 1.000 1.000 5.170 0.02411 *
## s(elev.mean):bb.spsoro 2.308 2.807 1.306 0.39447
## s(elev.mean):bb.spwurf 1.000 1.000 0.030 0.86283
## s(bb.transects) 1.000 1.000 3.695 0.05608 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.263 Deviance explained = 37.1%
## -REML = -11.722 Scale est. = 0.033787 n = 217print(plot(gam_d.2012, select = c(1, 3, 6, 8, 11)), pages = 1)
Let’s revisit option 1.
If we’re not committed to a story about elevation, can we explain our network metrics by including a larger suite of covariates, some of which may be collinear with elevation?
First, let’s inspect collinearity amongst the covariates. A PCA won’t be much help here since some of these relationships are nonlinear. So, we’ll make a pairs plot and interpret it manually. There’s a lot that could be said about these graphs, but here are some highlights.
- Elevation x covariates
- Roughly quadratic relationship with BB richness, BB abundance (log-transformed), and floral richness
- Roughly linear (negative) with floral abundance (log-transformed)
- More complex relationship with BB and floral transects
- Collinearity among covariates
- BB richness and BB abundance are very strongly correlated
- There is also a strong correlation between BB richness and floral richness
- In 2011, there is a very strong correlation between BB abundance and floral abundance; interestingly, this relationship does not hold up in the other two years
For a start, we will use all our covariates except for floral transects, since that one has to most tangential relationship to d’ (it’s only effect would be via the estimation of floral species’ relative abundances). We will use the select = TRUE option to perform double-penalty variable selection; this potentially penalizes some terms out of the model. 
Network/group-level metrics
For these models, we will use the bs = “fs” factor smooth method. This preserves degrees of freedom by forcing all smooths to have the same wiggliness.
Results
- Niche overlap (37% explained)
- Weakly quadratic negative effect of floral richness
- Linear positive effect of (log) BB abundance
- No significant effect of elevation
- Linear or weakly quadratic positive effect of transect count in 2011 and 2012, but negative (weakly quadratic) effect in 2010.
- C score (57% explained)
- Cubic negative effect of BB richness.
- Positive linear effect of elevation
- Negative linear effect of transect number
Discussion These models explain a lot more variance than the elevation-only models did. Interestingly, the effect of elevation drops out entirely in the niche overlap model with the addition of the other covariates. This suggests that the effects of elevation may, indeed, be mediated by effects on abundance and diversity, but moreover that abundance and diversity affect network metrics apart from the influence of elevation.
The story for niche overlap seems to be that diets become more similar as floral richness goes down and bumble bee abundance goes up. The relationship with floral richness seems easy to interpret: when there are fewer options, bumble bees are forced to overlap more in their choices. The positive relationship with BB abundance might be a samping effect; when a lot of bumble bees are seen, it is more likely that their observed diets will overlap.
C score increased lineary with elevation, indicating that, as one moves upslope, BB diets tend to become more disaggregated (in a presence/absence sense – C score is basically an inverse and binary version of niche overlap). It also decreased linearly with transect number, which is probably a sampling effect: with higher sampling intensity, more links are observed, resulting in more observed overlap and less disaggregation. The negative cubic relationship with BB richness – and remember that BB richness and BB abundance are almost perfectly correlated – also suggests a sampling effect: when there are more and more kinds of BB, observed overlap is higher.
# Niche overlap
gam_niche.overlap.global <- gam(niche.overlap.HL ~
s(bb.rich, k = 4) +
s(fl.rich, k = 4) +
s(log.fl.abund, k = 4) +
s(log.bb.abund, k= 4) +
s(elev.mean, year, bs = "fs", k = 5) +
s(bb.transects, year, bs = "fs", k = 4),
data = netmet,
family = "gaussian",
method = "REML",
select = TRUE) %>% getViz()
check.gamViz(gam_niche.overlap.global) ##
## Method: REML Optimizer: outer newton
## full convergence after 14 iterations.
## Gradient range [-1.034418e-05,3.131483e-05]
## (score -55.64349 & scale 0.009635063).
## Hessian positive definite, eigenvalue range [3.24396e-07,36.04731].
## Model rank = 40 / 40
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(bb.rich) 3.00e+00 1.27e-05 0.97 0.32
## s(fl.rich) 3.00e+00 2.11e+00 1.11 0.78
## s(log.fl.abund) 3.00e+00 6.61e-01 1.05 0.55
## s(log.bb.abund) 3.00e+00 8.84e-01 0.89 0.13
## s(elev.mean,year) 1.50e+01 1.29e+00 0.94 0.28
## s(bb.transects,year) 1.20e+01 1.67e+00 1.05 0.67
summary(gam_niche.overlap.global)##
## Family: gaussian
## Link function: identity
##
## Formula:
## niche.overlap.HL ~ s(bb.rich, k = 4) + s(fl.rich, k = 4) + s(log.fl.abund,
## k = 4) + s(log.bb.abund, k = 4) + s(elev.mean, year, bs = "fs",
## k = 5) + s(bb.transects, year, bs = "fs", k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.27462 0.01889 14.54 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(bb.rich) 1.268e-05 3 0.000 0.42032
## s(fl.rich) 2.113e+00 3 4.819 0.00059 ***
## s(log.fl.abund) 6.610e-01 3 0.650 0.06306 .
## s(log.bb.abund) 8.844e-01 3 2.550 0.00137 **
## s(elev.mean,year) 1.295e+00 14 0.158 0.13894
## s(bb.transects,year) 1.671e+00 10 0.456 0.03922 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.311 Deviance explained = 37.4%
## -REML = -55.643 Scale est. = 0.0096351 n = 73print(plot(gam_niche.overlap.global, allTerms = T), pages = 1)
plot(sm(gam_niche.overlap.global, 6)) +
l_fitLine(alpha = 0.6) 
# C score
gam_cscore.global <- gam(C.score.HL ~
s(bb.rich, k = 4) +
s(fl.rich, k = 4) +
s(log.fl.abund, k = 4) +
s(log.bb.abund, k= 4) +
s(elev.mean, year, bs = "fs", k = 5) +
s(bb.transects, year, bs = "fs", k = 4),
data = netmet,
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
check.gamViz(gam_cscore.global) ##
## Method: REML Optimizer: outer newton
## full convergence after 15 iterations.
## Gradient range [-5.128725e-06,1.300864e-05]
## (score -37.88911 & scale 0.01162369).
## Hessian positive definite, eigenvalue range [2.399802e-06,34.07957].
## Model rank = 40 / 40
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(bb.rich) 3.00 2.78 1.20 0.94
## s(fl.rich) 3.00 2.34 1.02 0.52
## s(log.fl.abund) 3.00 1.60 0.95 0.33
## s(log.bb.abund) 3.00 1.14 1.13 0.87
## s(elev.mean,year) 15.00 1.93 1.08 0.68
## s(bb.transects,year) 12.00 1.26 1.09 0.73
summary(gam_cscore.global)##
## Family: gaussian
## Link function: identity
##
## Formula:
## C.score.HL ~ s(bb.rich, k = 4) + s(fl.rich, k = 4) + s(log.fl.abund,
## k = 4) + s(log.bb.abund, k = 4) + s(elev.mean, year, bs = "fs",
## k = 5) + s(bb.transects, year, bs = "fs", k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4809 0.0178 27.02 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(bb.rich) 2.781 2.953 4.604 0.0108 *
## s(fl.rich) 2.337 2.707 3.205 0.0680 .
## s(log.fl.abund) 1.603 1.962 2.295 0.1399
## s(log.bb.abund) 1.138 1.246 0.321 0.5483
## s(elev.mean,year) 1.935 14.000 0.362 0.0457 *
## s(bb.transects,year) 1.256 10.000 0.397 0.0313 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.49 Deviance explained = 56.8%
## -REML = -37.889 Scale est. = 0.011624 n = 73print(plot(gam_cscore.global, allTerms = T), pages = 1)
plot(sm(gam_cscore.global, 4)) +
l_fitLine(alpha = 0.6)
plot(sm(gam_cscore.global, 5)) +
l_fitLine(alpha = 0.6) 
Species-level d’
Results (only significant terms shown)
- 2010 (41% explained)
- B. gerstaeckeri (positive, linear)
- B. hortorum (positive, linear)
- B. mendax (positive, linear)
- B. pratorum (positive, quadratic)
- B. psithyrus (positive, linear)
- B. pyrenaeus (negative, quadratic)
- B. soroensis (positive, quadratic)
- 2011 (58% explained)
- Negative effect of BB richness (weakly cubic)
- Positive effect of FL richness (cubic)
- B. sensu-strictu (positive, linear)
- B. hortorum (positive, linear)
- B. hypnorum (positive, linear)
- B. lapidarius (positive, linear)
- B. monticola (positive, logistic)
- B. pascuorum (positive, linear)
- B. pratorum (positive, logistic)
- B. psithyrus (positive, linear)
- B. pyrenaeus (positive, linear)
- B. soroensis (positive, linear)
- B. wurflenii (positive, linear)
- Transects (positive, quadratic)
- 2012 (43% explained)
- Positive effect of BB richness (linear)
- Negative effect of FL abundance (cubic)
- Negative effect of BB abundance (linear)
- B. hortorum (negative, linear)
- B. mendax (negative, linear)
- B. pascuorum (negative, linear)
- B/ pyrenaeus (negative, linear)
Discussion This is interesting. I do not fully understand why, but the addition of abundance and diversity covariates actually strengthens the apparent effects of elevation. It could also have something to do with setting select = TRUE. The change is most evident in 2011, where most BB species now show a significant relationship with elevation. The overall patterns are the same, though. Elevation has a positive effect on discrimination in 2010 and 2011, but this switches to a negative relationship in 2012.
Conclusion
### Species-specific models
gam_d.2010.global <- gam(d ~
s(bb.rich, k = 4) +
s(fl.rich, k = 4) +
s(log.fl.abund, k = 4) +
s(log.bb.abund, k= 4) +
bb.sp +
s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
s(bb.transects, k = 4),
data = filter(spmet, year == 2010 &
!bb.sp %in% c("humi", "lapi")),
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
gam_d.2010.global <- gam(d ~
s(bb.rich, k = 4) +
s(fl.rich, k = 4) +
s(log.fl.abund, k = 4) +
s(log.bb.abund, k= 4) +
s(elev.mean, bb.sp, bs = "fs", k = 5) +
s(bb.transects, k = 4),
data = filter(spmet, year == 2010 &
!bb.sp %in% c("humi", "lapi")),
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
check.gamViz(gam_d.2010.global) ##
## Method: REML Optimizer: outer newton
## full convergence after 17 iterations.
## Gradient range [-2.25218e-06,3.597569e-06]
## (score 20.06302 & scale 0.0569079).
## Hessian positive definite, eigenvalue range [9.160175e-07,91.72902].
## Model rank = 86 / 86
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(bb.rich) 3.00 1.00 1.04 0.65
## s(fl.rich) 3.00 1.00 1.05 0.74
## s(log.fl.abund) 3.00 2.54 1.07 0.78
## s(log.bb.abund) 3.00 1.00 1.08 0.88
## s(elev.mean,bb.sp) 70.00 8.92 1.08 0.82
## s(bb.transects) 3.00 1.48 1.06 0.79## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
summary(gam_d.2010.global)##
## Family: gaussian
## Link function: identity
##
## Formula:
## d ~ s(bb.rich, k = 4) + s(fl.rich, k = 4) + s(log.fl.abund, k = 4) +
## s(log.bb.abund, k = 4) + s(elev.mean, bb.sp, bs = "fs", k = 5) +
## s(bb.transects, k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.36559 0.02585 14.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(bb.rich) 1.000 1.000 0.019 0.89075
## s(fl.rich) 1.000 1.000 0.000 0.99621
## s(log.fl.abund) 2.541 2.852 1.542 0.16655
## s(log.bb.abund) 1.000 1.000 0.522 0.47099
## s(elev.mean,bb.sp) 8.918 66.000 0.285 0.00573 **
## s(bb.transects) 1.485 1.756 0.719 0.37097
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.118 Deviance explained = 19.3%
## -REML = 20.063 Scale est. = 0.056908 n = 189print(plot(gam_d.2010.global, allTerms = TRUE), pages = 1)
print(plot(gam_d.2010.global, select = c(6, 7, 10, 14, 15, 16, 17)), pages = 1)
gam_d.2011.global <- gam(d ~
s(bb.rich, k = 4) +
s(fl.rich, k = 4) +
s(log.fl.abund, k = 4) +
s(log.bb.abund, k= 4) +
bb.sp +
s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
s(bb.transects, k = 4),
data = filter(spmet, year == 2011 &
!bb.sp %in% c("humi", "gers", "muci")),
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
gam_d.2011.global <- gam(d ~
s(bb.rich, k = 4) +
s(fl.rich, k = 4) +
s(log.fl.abund, k = 4) +
s(log.bb.abund, k= 4) +
s(elev.mean, bb.sp, bs = "fs", k = 5) +
s(bb.transects, k = 4),
data = filter(spmet, year == 2011 &
!bb.sp %in% c("humi", "gers", "muci")),
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
check.gamViz(gam_d.2011.global) ##
## Method: REML Optimizer: outer newton
## full convergence after 17 iterations.
## Gradient range [-5.129729e-05,0.0001512141]
## (score -52.69027 & scale 0.02657087).
## Hessian positive definite, eigenvalue range [2.727704e-06,114.903].
## Model rank = 81 / 81
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(bb.rich) 3.00 1.00 0.79 <2e-16 ***
## s(fl.rich) 3.00 1.00 0.79 <2e-16 ***
## s(log.fl.abund) 3.00 2.20 0.79 <2e-16 ***
## s(log.bb.abund) 3.00 1.66 0.77 <2e-16 ***
## s(elev.mean,bb.sp) 65.00 18.68 0.79 <2e-16 ***
## s(bb.transects) 3.00 2.50 0.81 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
summary(gam_d.2011.global)##
## Family: gaussian
## Link function: identity
##
## Formula:
## d ~ s(bb.rich, k = 4) + s(fl.rich, k = 4) + s(log.fl.abund, k = 4) +
## s(log.bb.abund, k = 4) + s(elev.mean, bb.sp, bs = "fs", k = 5) +
## s(bb.transects, k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3782 0.0297 12.73 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(bb.rich) 1.000 1.000 7.494 0.00672 **
## s(fl.rich) 1.000 1.000 1.204 0.27386
## s(log.fl.abund) 2.204 2.612 2.502 0.12074
## s(log.bb.abund) 1.663 2.014 1.541 0.22276
## s(elev.mean,bb.sp) 18.682 64.000 1.589 3.91e-13 ***
## s(bb.transects) 2.501 2.810 5.521 0.00660 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.421 Deviance explained = 48.8%
## -REML = -52.69 Scale est. = 0.026571 n = 235print(plot(gam_d.2011.global, allTerms = TRUE), pages = 1)
print(plot(gam_d.2011.global, select = c(1, 2, 5, 6, 7, 9, 11, 12, 13, 14, 15, 16, 17, 18)), pages = 1)
gam_d.2012.global <- gam(d ~
s(bb.rich, k = 4) +
s(fl.rich, k = 4) +
s(log.fl.abund, k = 4) +
s(log.bb.abund, k= 4) +
bb.sp +
s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
s(bb.transects, k = 4),
data = filter(spmet, year == 2012 &
!bb.sp %in% c("humi", "hypn", "muci")),
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
gam_d.2012.global <- gam(d ~
s(bb.rich, k = 4) +
s(fl.rich, k = 4) +
s(log.fl.abund, k = 4) +
s(log.bb.abund, k= 4) +
s(elev.mean, bb.sp, bs = "fs", k = 5) +
s(bb.transects, k = 4),
data = filter(spmet, year == 2012 &
!bb.sp %in% c("humi", "hypn", "muci")),
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
gam_d.2012.global <- gam(d ~
s(bb.rich, k = 4) +
s(fl.rich, k = 4) +
s(log.fl.abund, k = 4) +
s(log.bb.abund, k= 4) +
s(elev.mean, bb.sp, bs = "fs", k = 5) +
s(bb.transects, k = 4),
data = spmet,
family = "gaussian",
method = "REML",
select = FALSE) %>% getViz()
check.gamViz(gam_d.2012.global) ##
## Method: REML Optimizer: outer newton
## full convergence after 17 iterations.
## Gradient range [-1.927926e-05,9.303341e-05]
## (score -45.08819 & scale 0.04610161).
## Hessian positive definite, eigenvalue range [1.798701e-06,325.6101].
## Model rank = 96 / 96
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(bb.rich) 3.00 1.00 0.89 <2e-16 ***
## s(fl.rich) 3.00 2.11 0.89 <2e-16 ***
## s(log.fl.abund) 3.00 1.35 0.89 <2e-16 ***
## s(log.bb.abund) 3.00 1.00 0.88 0.005 **
## s(elev.mean,bb.sp) 80.00 11.84 0.87 <2e-16 ***
## s(bb.transects) 3.00 1.00 0.88 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
summary(gam_d.2012.global)##
## Family: gaussian
## Link function: identity
##
## Formula:
## d ~ s(bb.rich, k = 4) + s(fl.rich, k = 4) + s(log.fl.abund, k = 4) +
## s(log.bb.abund, k = 4) + s(elev.mean, bb.sp, bs = "fs", k = 5) +
## s(bb.transects, k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.37188 0.02279 16.32 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(bb.rich) 1.000 1.000 0.021 0.8844
## s(fl.rich) 2.114 2.522 2.536 0.0724 .
## s(log.fl.abund) 1.346 1.612 6.837 0.0121 *
## s(log.bb.abund) 1.000 1.000 0.305 0.5808
## s(elev.mean,bb.sp) 11.841 77.000 0.870 1.21e-10 ***
## s(bb.transects) 1.000 1.000 3.882 0.0492 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.123 Deviance explained = 14.7%
## -REML = -45.088 Scale est. = 0.046102 n = 657print(plot(gam_d.2012.global, select = c(1, 3, 4, 6, 10, 12, 15)), pages = 1)
plot(gam_d.2012.global, select = c(5))
plot(gam_d.2011.global, select = c(5))
plot(gam_d.2010.global, select = c(5))
Conclusion
When niche overlap and C score are modeled only as functions of elevation and transect count, we can find interpretable, if weak, elevational patterns. When we include the other covariates, though, the influence of elevation is largely eclipsed by that of BB abundance/richness, and I am not confident that these relationships can be interpreted as anything more than sampling effects. This casts doubt on the legitimacy of the elevation patterns found when abundance and richness covariates are excluded; it’s hard to make the case that they are not just sampling effects, too.
The preliminary analysis below are not finished yet, so I won’t include the results at this stage.
Module analysis
First for the metaweb